/*
 * Utilities.cpp
 *
 *  Created on: 10 Aug 2011
 *      Author: Allan
 */

#include "Utilities.h"

// C++ includes
#include <iostream>

// Boost includes
#include <boost/foreach.hpp>

// Eigen includes
#include <Eigen/Core>

const vector<EquilibriumCondition> ConvertEquilibriumConditions(const vector<string>& equilibriumConditions, const Multiphase& multiphase, const ReactionSystem& equilibriumReactions)
{
	vector<EquilibriumCondition> eConditions;
	
	BOOST_FOREACH(const string& condition, equilibriumConditions)
	{
		EquilibriumCondition eCondition(multiphase, equilibriumReactions);
	
		vector<string> words = SplitLineAndClean(condition, ":= ");
	
		if(words[0] == "ChargeBalanced")
			eCondition.SetChargeBalance();
	
		else
		if(words[0] == "Composition")
		{
			const string species = words[1]; 
			const double moles = atof(words[2].c_str()); 
			
			eCondition.SetComposition(species, moles);
		}
	
		else
		if(words[0] == "TotalComposition")
		{
			const string species = words[1];
			const double moles = atof(words[2].c_str());
		
			eCondition.SetTotalComposition(species, moles);
		}
	
		else
		if(words[0] == "Activity")
		{
			const string species = words[1];
			const double activity = atof(words[2].c_str());
	
			eCondition.SetActivity(species, activity);
		}
	
		else
		if(words[0] == "pH")
		{
			const double pH = atof(words[1].c_str());
	
			eCondition.SetAcidity(pH);
		}
	
		else
		if(words[0] == "PartialPressure")
		{
			const string gas = words[1];
			const double Pg = atof(words[2].c_str());
	
			eCondition.SetPartialPressure(gas, Pg);
		}
	
		else cerr << "Error: The specified equilibrium condition " << condition << " is incorrect or invalid." << endl;
		
		eConditions.push_back(eCondition);
	}
	
	return eConditions;
}

const VectorXd GuessSpeciesComposition(const vector<string>& generalConditions, const Multiphase& multiphase, double T, double P)
{
	VectorXd n = VectorXd::Constant(multiphase.GetSpecies().size(), 1.0E-7);
	
	const SpeciesMap& speciesMap = multiphase.GetSpeciesMap();
	
	BOOST_FOREACH(const string& condition, generalConditions)
	{
		vector<string> words = SplitLineAndClean(condition, ":= ");
	
		if(words[0] == "Composition")
		{
			const string species = words[1]; 
			const double moles = atof(words[2].c_str()); 
			
			n[speciesMap[species]] = moles;
		}
	
		if(words[0] == "TotalComposition")
		{
			const string species = words[1];
			const double moles = atof(words[2].c_str());
		
			n[speciesMap[species]] = (moles > 0) ? moles : 1.0E-7;
		}
	
		if(words[0] == "Activity")
		{
			const string species = words[1];
			const double activity = atof(words[2].c_str());
	
			n[speciesMap[species]] = activity;
		}
	
		if(words[0] == "pH")
		{
			const double pH = atof(words[1].c_str());
	
			n[speciesMap["H[+]"]] = std::pow(10.0, -pH);
		}
	
		if(words[0] == "PartialPressure")
		{
			const string gas = words[1];
			const double Pg = atof(words[2].c_str());
	
			n[speciesMap[gas]] = Pg/P;
		}
	}
	
	return n;
}

const vector<string> ExtractEquilibriumConditions(const vector<string>& generalConditions, const Multiphase& multiphase, const ReactionSystem& equilibriumReactions)
{
	vector<string> equilibriumConditions;
	
	const vector<string> jSpecies = equilibriumReactions.Products();
	const vector<string> eSpecies = equilibriumReactions.Reactants();
	
	auto isPrimary     = [&](const string& species) { return std::find(jSpecies.begin(), jSpecies.end(), species) != jSpecies.end(); };
	auto isEquilibrium = [&](const string& species) { return std::find(eSpecies.begin(), eSpecies.end(), species) != eSpecies.end(); };
	
	BOOST_FOREACH(const string& condition, generalConditions)
	{
		vector<string> words = SplitLineAndClean(condition, ":= ");
	
		if(words[0] == "Composition" || words[0] == "TotalComposition")
		{
			const string species = words[1];
			
			if(isPrimary(species) || isEquilibrium(species))
				equilibriumConditions.push_back(condition);
		}
		else
			equilibriumConditions.push_back(condition);
	}
	
	return equilibriumConditions;
}
